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A fundamental macroscopic description of a magnetized plasma is the Vlasov equation supplemented by 
the nonlinear inverse-square force Fokker-Planck collision operator [Rosenbluth et al., Phys. Rev., 107 , 1957], 
The Vlasov part describes advection in a six-dimensional phase space whereas the collision operator involves 
friction and diffusion coefficients that are weighted velocity-space integrals of the particle distribution function. 
The Fokker-Planck collision operator is an integro-differential, bilinear operator, and numerical discretization 
of the operator is far from trivial. In this letter, we describe a new approach to discretize the entire kinetic 
system based on an expansion in Gaussian Radial Basis functions (RBFs). This approach is particularly well- 
suited to treat the collision operator because the friction and diffusion coefficients can be analytically calculated. 
Although the RBF method is known to be a powerful scheme for the interpolation of scattered multidimensional 
data, Gaussian RBFs also have a deep physical interpretation as local thermodynamic equilibria. In this letter we 
outline the general theory, highlight the connection to plasma fluid theories, and also give 2D and 3D numerical 
solutions of the nonlinear Fokker-Planck equation. A broad spectrum of applications for the new method is 
anticipated in both astrophysical and laboratory plasmas. 


Motivation A fundamental macroscopic description of a 
magnetized plasma is the Vlasov equation supplemented by 
the nonlinear inverse-square force Fokker-planck collision op¬ 
erator [1] 
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where f a is the distribution of species a with charge e a and 
mass m a . The electric and magnetic fields depend on the 
dsitribution f„ through Maxwell’s equations. This model as¬ 
sumes a statistical description of Coulomb interaction in the 
limit of small-angle scattering, with the changes in f a due to 
collisions with species b described by 
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The friction and diffusion coefficients are given by the expres¬ 
sions 
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where L a b = {e a eb/rn a £o) 2 In A a b and lnA ai) is the Coulomb 
logarithm which respresents a physical cut-off for small-angle 
collisions. The Rosenbluth potentials appearing in the friction 
and diffusion coefficients are weighted integrals of the distri¬ 
bution function 
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and they satisfy the velocity-space Poisson equations 'V 2 tpb = 
tp b and Vyipb = fb • Expressed in terms of the potential func¬ 
tions, the Fokker-Planck collision operator is 
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where p ab = m a /m b — 1. A common approach for numeri¬ 
cal evaluation of C ab follows a two-phase method where one 
first inverts the velocity-space Laplacian operators and then 
directly evaluates the collision operator. Boundary conditions 
for the Poisson equations can be obtained by limiting the so¬ 
lution to a sub-space and evaluating the expressions at the 
boundary using a multipole expansion of the potentials [2], or 
by imposing the free-space solution outside the sub-space [3], 
Another sophisticated approach is based on fast spectral de¬ 
composition as described in Ref. [4], In this letter, we propose 
a new approach using a mesh-free shifted-Maxwellian repre¬ 
sentation which is intuitively appealing and straightforward 
to implement. The solution thus obtained is C°° smooth, ex¬ 
tends to v —t 00, and allows compact representation of any in¬ 
teresting macroscopic quantity (number, momentum, energy 
density, and so on). 

The Gaussian RBF Method To solve the kinetic equation, 
Eq. (1), we write the total distribution as a finite sum of shifted 
Maxwellians 

fa (x, V, t) = W l a (x, t) Fl (x, v) , 

where each Maxwellian, F* = (ya/n) 3 ^ 2 exp [—7* (v — vl) 2 ], 
is normalized to unity and the weights vj'„ allowed to evolve in 
time. The width parameters 7* and mean velocities v* can be 
arbitrary functions of position. The shifted Maxwellians are 
nothing other than Gaussian Radial Basis functions (RBFs) 
which have found numerous applications in applied mathe¬ 
matics - in particular for the construction of neural networks 
[5]. For compactness, in what follows we will retain the spa¬ 
tial dependence of all quantities but will not write the depen¬ 
dence explicitly. The potential functions then take the form 
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in terms of the functions <f>(s) = erf(s)/s and 'k(s) = [s + 
l/(2s)]erf(s) + exp(— s 2 )/ yTr, where erf(s) is the error func¬ 
tion. We thus find a simple bilinear expression for the com¬ 
plete nonlinear operator 

Cab = £ (f ) (t) Cab 0) , 
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where the Gaussian RBF collision tensor is 
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such that C kk (v) = 0. As we have analytical expressions for 
Fa{\), y>l(v), and t/j“(v), the tensor C k l(y) is easy to imple¬ 
ment and fast to evaluate at any point in velocity space. In 
problems with azimuthal symmetry, a 2D RBF scheme can be 
developed with axisymmetric ring -like RBF-basis: 

Fi = ( 7 i/ 7 r) 3 / 2 / 0 (2 7i « i ,x«x)e- 7i(,J l|-^lP^^i+^A) , 


where I 0 is the order-zero modified Bessel function of the first 
kind and (iJ_ L ,tj|j) are the cylindrical velocity-space coordi¬ 
nates. Although explicit expressions for axisymmetric RBF 
potentials <pi and ipi are not available in a closed form, they 
can be evaluated numerically to machine precision at any re¬ 
quested point. 

Collocation Options We describe two different methods 
for obtaining an ordinary differential equation for the time 
evolution of weights: the Galerkin and the center-collocation 
projections. In the Galerkin method, the kinetic equation - 
already expanded in the RBF basis - is multiplied by each 
individual basis function and then integrated over the entire 
domain. Conversely, in the center-collocation method, the ki¬ 
netic equation is evaluated at the center of each RBF. Both 
methods yield N equations for the N RBF weights, w l a , and 
result in a differential equation for the weights that can be 
written in a matrix form. For the moment, let us illustrate the 
method by considering the spatially homogeneous case with 
no electromagnetic fields. Then, the matrix equation is 


for the Galerkin and center-collocation, respectively. Obtain¬ 
ing the center-collocation tensor (C*^)cc is merely a matter 
of evaluating the RBF tensor C k i(y) at the collocation points. 
Evaluation of the tensor (C*^)gp for the Galerkin projection 
is somewhat more complicated, although the result may be 
potentially be more accurate or robust. Nevertheless, to main¬ 
tain simplicity, we focus hereafter on the center collocation- 
method and omit the CC subscript for brevity. In this case the 
RBF equations for the full nonlinear system become 

^2 Mij Lijwi = £ w k a wic l ab , V i € 1, 2, 3,..., (3) 
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where the operator £ 

Cij “ §~ t + Vi ' V + ^ [(Vj - Vi) ■ E + (vj x Vi) • b] 

retains the familiar appearance of the Vlasov operator even 
though the velocity-space has been completely removed from 
the problem. Note that £ depends explicitly on species a and 
implicitly on b via the Maxwell equations. In £ we have de¬ 
fined the temperature-like parameter of = m a /(2 , yi) and also 
dropped some additional terms that arise if the parameters 7 i 
and Vi depend on position. The RBF-kinetic equation (3) de¬ 
scribes collisional fluid-like evolution of the weights in time 
and space. One physically appealing feature of the RBF ex¬ 
pansion is that familiar expressions for macroscopic fluid mo¬ 
ments retain their intuitive form: 
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In the Galerkin projection (GP), the matrix M,, is symmetric, 
typically diagonally dominant, and given by the expression 
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whereas in the center-collocation (CC) method, the matrix 
Mij is no longer necessarily symmetric, but can still be domi¬ 
nated by the diagonal components: 

(M^)cc = fj(\i ) = (7 j/A 312 exp[— 7 j(v; - \j) 2 } . 


On the right-hand-side of Eq. (2), the tensor C a k h f becomes 


id “)gp = J dYfi(Y) Cab (v) , (Clf)CC = C££(v<) , 


Another is that Maxwell’s equations can also written com¬ 
pactly: 

V • E = 47 t £ e a n a = 4tt £ e a w l a , 
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where for simplicity we have neglected the displacement cur¬ 
rent. In contrast to the standard fluid approach, where calcula¬ 
tion of the higher-order velocity-space moments becomes in¬ 
creasingly cumbersome, the RBF-kinetic approach, however, 
offers a straightforward and intuitive alternative by describing 
the velocity space with Gaussian functions that naturally ap¬ 
pear in the kinetic theory as thermodynamic equilibrium solu¬ 
tions. Moreover, so long as the RBF spacing is equal to about 
one e-folding length, the problem remains well-conditioned 
even for hundreds or thousands of basis functions. 
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FIG. 1. Convergence of a steady-state collision operator. The lines 
with square (star) markers represent the 3D (axisymmetric 2D) solu¬ 
tions and A is the distance between adjacent collocation points. 


FIG. 2. Convergence of the velocity space moments in a relaxation 
study. Here 7 = 0.5/A 2 was used and A is the distance between 
adjacent collocation points. 


It should also be evident that the Gaussian RBF method is 
well-suited to linearized models.The linearized collision oper¬ 
ator in this case is just the sum of C° l and C l °; namely, the first 
row and column of the nonlinear collision tensor. Because the 
method is mesh-free, it is also natually suited to multi-species 
plasmas, unlike contemporary algorithms that are based on a 
velocity mesh [ 6 ], 

Nonlinear Simulations in 2D and 3D To assess the 
viability of the Gaussian RBF approach, we study a single 
species plasma using a uniform grid for the collocation points 
v, and a fixed value 7 ,; = 7 for all RBFs. We also normal¬ 
ize time with the collision time, i.e., r = Our problem is 
thus to solve the nonlinear single species Fokker-Planck equa¬ 
tion 


d£ = ,2 _ 9 2 tp _ d 2 f = 
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We emphasize that this is a highly nontrivial numerical prob¬ 
lem. First, we test how well the discretization preserves a 
(stationary) Maxwellian equilibrium state, where the temper¬ 
ature of the equilibrium is fixed and much warmer than the 
individual RBFs. For the analytic equilibrium we choose 
f M = (/3/7t) 3/2 exp [—,0v 2 ], with (5 = 1/4 and arrange the 
collocation points into uniform rectilinear grids in spaces 
[v x , v y , v z ] £ [—6.5,6.5] x [—6.5, 6.5] x [—6.5, 6.5] and [v||,'Ux] £ 
[—6.5,6.5] x [0,6.5] for the full 3D and axisymmetric 2D meth¬ 
ods respectively. The equilibrium is projected onto the RBF 
basis to obtain the numerical equivalent of the steady state 
distribution function. Then the collision operator in Eq. (4) 
is evaluated for an increasing number of collocation points 
and different values of 7 . As an error metric, we choose 
ma x{C[f]/f}, presented in Fig. 1. This test is sensitive to 
the error both in the region interior to the RBF collocation 
centers, as well as to the region beyond them. As is evident 
from Fig. 1, the maximum of the global value decreases when 
the number of collocation points is increased indicating that 
a numerical equivalent to the analytical steady-state can be 
reached. 


Next, we follow the relaxation of a non-trivial distribu¬ 
tion function to the equilibrium state and track the conser¬ 
vation of number, momentum, and energy densities. The 
initial state of the distribution function is set to /(v, r = 
0) = X)Li e -/3i{v-Vi)2 , where ft * = 1/5, Vi = (3,0,0), 
V 2 = (—3,0,0) for the ( v x ,v y ,v z ) components and respec¬ 
tively. The initial state is thus axisymmetric with respect to 
the axis v x and the the axisymmetric setting is thus chosen to 
be (17 = v x , v\ = v 2 + v 2 ). The collocation points are chosen 
uniformly in the regions [v x , v y ,v z \ £ [— 8 , 8 ] x [— 8 , 8 ] x [— 8 , 8 ] 
and [w||, v_l] £ [— 12 , 12 ] x [ 0 , 12 ], the initial state projected to 
the Gaussian RBF basis to get initial values for the weights, 
and the weights then evolved in time according to the Eq. (2) 
using a standard fourth-order Runge-Kutta method. From the 
time-evolution of the weights, we have calculated the evolu¬ 
tion of the velocity-space moments, and recorded their max¬ 
imal deviations from the initial values during the simulation. 
The results are illustrated in Fig. 2, where one observes good 
conservation of number (n), momentum (velocity components 
v x , v y , and v z ), and energy (£,) densities, and convergence of 
their relative error towards zero when the number of RBFs is 
increased for both full 3D and axisymmetric 2D solvers. Even 
better conservation properties can be reached by embedding 
the conservation laws into the time-evolution of the weights. 
For example, replacing one row in the matrix (M^ jcc with 
ones (unity moments of the basis functions) and setting the 
corresponding element in Ciki,cc to zero neglects information 
from one of the basis nodes but adds the density conservation. 
With this approach we observed, e.g., density conservation up 
to 12 digits without relevant increase in the computation time. 

The time evolution of the distribution in the conservation 
study using 15 3 3D RBFs and 30 x 15 axisymmetric RBFs is 
illustrated in Fig. 3 with slices in (v x , ly/)-plane for the full 3D 
solver and in ( 17 , i>x (-plane for the 2D axisymmetric solver. 
The time slices are chosen for the initial state (r = 0) in 
Figs. 3(a) and 3(e), for the beginning of the relaxation process 
(r = 3) in Figs. 3(b) and 3(f), for the merging phase towards 
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FIG. 3. Relaxation of an axisymmetric initial state towards spherically symmetric equilibrium state from the 3D solver (a,b,c,d) and axisym- 
metric 2D solver (e,f,g,h). The solution is calculated using 15 x 15 x 15 3D RBFs and 30 x 15 axisymmetric RBFs. 



FIG. 4. Time evolution of the u 4 -moment during the relaxation study 
(the almost overlapping red and green curves). The lower (upper) 
dashed line represents the analytical initial (equilibrium) values. 


the equilibrium (r = 6) in Figs. 3(c) and 3(g), and for the final 
state (r = 40) where the distribution function is close to spher¬ 
ical symmetric and equilibrium in Figs. 3(d) and 3(h). Both 
the 2D and 3D solvers qualitatively describe the same solution 
to the initially axisymmetric problem and preserve the initial 
symmetry. As a more quantitative measure, we also follow 
the time evolution of the //'-moment of the distribution func¬ 
tion which is not a conserved quantity. From the analytical 
initial state one can calculate the value to be 4.98 x 10 4 and 
for the corresponding analytical equilibrium state the value is 
5.65 x 10 4 . The time evolution of v 4 -moment from the numer¬ 
ical calculation is shown in Fig. 4. As is clear, the 2D and 


3D solutions agree with each other and also confirm our claim 
that the solution in Figs. 3(d) and 3(h) is close to the analytical 
equilibrium state. In other words, we have solved the relax¬ 
ation problem accurately in both 2D and 3D and demonstrated 
the capabilities of the RBFs for discretizing the nonlinear col¬ 
lision operator. 

Summary In this letter, we have derived a completely new 
approach to collisional dynamics in plasmas. We have shown 
how to implement a Gaussian radial-basis-function ansatz to 
discretize the velocity space, and demonstrated its capabili¬ 
ties in non-linear Fokker-Planck calculations in both full 3D 
and axisymmetric 2D velocity-space. In a wider connection to 
kinetic theory we have provided also the discretization of the 
advection terms and given analytical ansatz expressions for all 
relevant fluid-moments. In future, our method could be used 
to solve various problems where non-equilibrium kinetic ef¬ 
fects are important, as long as the Landau-approximation for 
the collision integral is valid. 
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